<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>sstebz.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SSTEBZ.1"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a>( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
     $                   M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">     8-18-00:  Increase FUDGE factor for T3E (eca)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          ORDER, RANGE
      INTEGER            IL, INFO, IU, M, N, NSPLIT
      REAL               ABSTOL, VL, VU
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * )
      REAL               D( * ), E( * ), W( * ), WORK( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="SSTEBZ.23"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a> computes the eigenvalues of a symmetric tridiagonal
</span><span class="comment">*</span><span class="comment">  matrix T.  The user may ask for all eigenvalues, all eigenvalues
</span><span class="comment">*</span><span class="comment">  in the half-open interval (VL, VU], or the IL-th through IU-th
</span><span class="comment">*</span><span class="comment">  eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  To avoid overflow, the matrix must be scaled so that its
</span><span class="comment">*</span><span class="comment">  largest element is no greater than overflow**(1/2) *
</span><span class="comment">*</span><span class="comment">  underflow**(1/4) in absolute value, and for greatest
</span><span class="comment">*</span><span class="comment">  accuracy, it should not be much smaller than that.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  See W. Kahan &quot;Accurate Eigenvalues of a Symmetric Tridiagonal
</span><span class="comment">*</span><span class="comment">  Matrix&quot;, Report CS41, Computer Science Dept., Stanford
</span><span class="comment">*</span><span class="comment">  University, July 21, 1966.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RANGE   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'A': (&quot;All&quot;)   all eigenvalues will be found.
</span><span class="comment">*</span><span class="comment">          = 'V': (&quot;Value&quot;) all eigenvalues in the half-open interval
</span><span class="comment">*</span><span class="comment">                           (VL, VU] will be found.
</span><span class="comment">*</span><span class="comment">          = 'I': (&quot;Index&quot;) the IL-th through IU-th eigenvalues (of the
</span><span class="comment">*</span><span class="comment">                           entire matrix) will be found.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ORDER   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'B': (&quot;By Block&quot;) the eigenvalues will be grouped by
</span><span class="comment">*</span><span class="comment">                              split-off block (see IBLOCK, ISPLIT) and
</span><span class="comment">*</span><span class="comment">                              ordered from smallest to largest within
</span><span class="comment">*</span><span class="comment">                              the block.
</span><span class="comment">*</span><span class="comment">          = 'E': (&quot;Entire matrix&quot;)
</span><span class="comment">*</span><span class="comment">                              the eigenvalues for the entire matrix
</span><span class="comment">*</span><span class="comment">                              will be ordered from smallest to
</span><span class="comment">*</span><span class="comment">                              largest.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the tridiagonal matrix T.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VL      (input) REAL
</span><span class="comment">*</span><span class="comment">  VU      (input) REAL
</span><span class="comment">*</span><span class="comment">          If RANGE='V', the lower and upper bounds of the interval to
</span><span class="comment">*</span><span class="comment">          be searched for eigenvalues.  Eigenvalues less than or equal
</span><span class="comment">*</span><span class="comment">          to VL, or greater than VU, will not be returned.  VL &lt; VU.
</span><span class="comment">*</span><span class="comment">          Not referenced if RANGE = 'A' or 'I'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IL      (input) INTEGER
</span><span class="comment">*</span><span class="comment">  IU      (input) INTEGER
</span><span class="comment">*</span><span class="comment">          If RANGE='I', the indices (in ascending order) of the
</span><span class="comment">*</span><span class="comment">          smallest and largest eigenvalues to be returned.
</span><span class="comment">*</span><span class="comment">          1 &lt;= IL &lt;= IU &lt;= N, if N &gt; 0; IL = 1 and IU = 0 if N = 0.
</span><span class="comment">*</span><span class="comment">          Not referenced if RANGE = 'A' or 'V'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ABSTOL  (input) REAL
</span><span class="comment">*</span><span class="comment">          The absolute tolerance for the eigenvalues.  An eigenvalue
</span><span class="comment">*</span><span class="comment">          (or cluster) is considered to be located if it has been
</span><span class="comment">*</span><span class="comment">          determined to lie in an interval whose width is ABSTOL or
</span><span class="comment">*</span><span class="comment">          less.  If ABSTOL is less than or equal to zero, then ULP*|T|
</span><span class="comment">*</span><span class="comment">          will be used, where |T| means the 1-norm of T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          Eigenvalues will be computed most accurately when ABSTOL is
</span><span class="comment">*</span><span class="comment">          set to twice the underflow threshold 2*<a name="SLAMCH.82"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>('S'), not zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The n diagonal elements of the tridiagonal matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E       (input) REAL array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          The (n-1) off-diagonal elements of the tridiagonal matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  M       (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The actual number of eigenvalues found. 0 &lt;= M &lt;= N.
</span><span class="comment">*</span><span class="comment">          (See also the description of INFO=2,3.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NSPLIT  (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of diagonal blocks in the matrix T.
</span><span class="comment">*</span><span class="comment">          1 &lt;= NSPLIT &lt;= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  W       (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On exit, the first M elements of W will contain the
</span><span class="comment">*</span><span class="comment">          eigenvalues.  (<a name="SSTEBZ.100"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a> may use the remaining N-M elements as
</span><span class="comment">*</span><span class="comment">          workspace.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IBLOCK  (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          At each row/column j where E(j) is zero or small, the
</span><span class="comment">*</span><span class="comment">          matrix T is considered to split into a block diagonal
</span><span class="comment">*</span><span class="comment">          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
</span><span class="comment">*</span><span class="comment">          block (from 1 to the number of blocks) the eigenvalue W(i)
</span><span class="comment">*</span><span class="comment">          belongs.  (<a name="SSTEBZ.108"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a> may use the remaining N-M elements as
</span><span class="comment">*</span><span class="comment">          workspace.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ISPLIT  (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The splitting points, at which T breaks up into submatrices.
</span><span class="comment">*</span><span class="comment">          The first submatrix consists of rows/columns 1 to ISPLIT(1),
</span><span class="comment">*</span><span class="comment">          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
</span><span class="comment">*</span><span class="comment">          etc., and the NSPLIT-th consists of rows/columns
</span><span class="comment">*</span><span class="comment">          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
</span><span class="comment">*</span><span class="comment">          (Only the first NSPLIT elements will actually be used, but
</span><span class="comment">*</span><span class="comment">          since the user cannot know a priori what value NSPLIT will
</span><span class="comment">*</span><span class="comment">          have, N words must be reserved for ISPLIT.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) REAL array, dimension (4*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace) INTEGER array, dimension (3*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  some or all of the eigenvalues failed to converge or
</span><span class="comment">*</span><span class="comment">                were not computed:
</span><span class="comment">*</span><span class="comment">                =1 or 3: Bisection failed to converge for some
</span><span class="comment">*</span><span class="comment">                        eigenvalues; these eigenvalues are flagged by a
</span><span class="comment">*</span><span class="comment">                        negative block number.  The effect is that the
</span><span class="comment">*</span><span class="comment">                        eigenvalues may not be as accurate as the
</span><span class="comment">*</span><span class="comment">                        absolute and relative tolerances.  This is
</span><span class="comment">*</span><span class="comment">                        generally caused by unexpectedly inaccurate
</span><span class="comment">*</span><span class="comment">                        arithmetic.
</span><span class="comment">*</span><span class="comment">                =2 or 3: RANGE='I' only: Not all of the eigenvalues
</span><span class="comment">*</span><span class="comment">                        IL:IU were found.
</span><span class="comment">*</span><span class="comment">                        Effect: M &lt; IU+1-IL
</span><span class="comment">*</span><span class="comment">                        Cause:  non-monotonic arithmetic, causing the
</span><span class="comment">*</span><span class="comment">                                Sturm sequence to be non-monotonic.
</span><span class="comment">*</span><span class="comment">                        Cure:   recalculate, using RANGE='A', and pick
</span><span class="comment">*</span><span class="comment">                                out eigenvalues IL:IU.  In some cases,
</span><span class="comment">*</span><span class="comment">                                increasing the PARAMETER &quot;FUDGE&quot; may
</span><span class="comment">*</span><span class="comment">                                make things work.
</span><span class="comment">*</span><span class="comment">                = 4:    RANGE='I', and the Gershgorin interval
</span><span class="comment">*</span><span class="comment">                        initially used was too small.  No eigenvalues
</span><span class="comment">*</span><span class="comment">                        were computed.
</span><span class="comment">*</span><span class="comment">                        Probable cause: your machine has sloppy
</span><span class="comment">*</span><span class="comment">                                        floating-point arithmetic.
</span><span class="comment">*</span><span class="comment">                        Cure: Increase the PARAMETER &quot;FUDGE&quot;,
</span><span class="comment">*</span><span class="comment">                              recompile, and try again.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Internal Parameters
</span><span class="comment">*</span><span class="comment">  ===================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RELFAC  REAL, default = 2.0e0
</span><span class="comment">*</span><span class="comment">          The relative tolerance.  An interval (a,b] lies within
</span><span class="comment">*</span><span class="comment">          &quot;relative tolerance&quot; if  b-a &lt; RELFAC*ulp*max(|a|,|b|),
</span><span class="comment">*</span><span class="comment">          where &quot;ulp&quot; is the machine precision (distance from 1 to
</span><span class="comment">*</span><span class="comment">          the next larger floating point number.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  FUDGE   REAL, default = 2
</span><span class="comment">*</span><span class="comment">          A &quot;fudge factor&quot; to widen the Gershgorin intervals.  Ideally,
</span><span class="comment">*</span><span class="comment">          a value of 1 should work, but on machines with sloppy
</span><span class="comment">*</span><span class="comment">          arithmetic, this needs to be larger.  The default for
</span><span class="comment">*</span><span class="comment">          publicly released versions should be large enough to handle
</span><span class="comment">*</span><span class="comment">          the worst machine around.  Note that this has no effect
</span><span class="comment">*</span><span class="comment">          on accuracy of the solution.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE, TWO, HALF
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
     $                   HALF = 1.0E0 / TWO )
      REAL               FUDGE, RELFAC
      PARAMETER          ( FUDGE = 2.1E0, RELFAC = 2.0E0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            NCNVRG, TOOFEW
      INTEGER            IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
     $                   IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
     $                   ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL,
     $                   NWU
      REAL               ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN,
     $                   TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            IDUMMA( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.193"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.194"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
      REAL               <a name="SLAMCH.195"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="LSAME.196"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.196"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>, <a name="SLAMCH.196"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="SLAEBZ.199"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>, <a name="XERBLA.199"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode RANGE
</span><span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.210"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'A'</span> ) ) THEN
         IRANGE = 1
      ELSE IF( <a name="LSAME.212"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'V'</span> ) ) THEN
         IRANGE = 2
      ELSE IF( <a name="LSAME.214"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'I'</span> ) ) THEN
         IRANGE = 3
      ELSE
         IRANGE = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode ORDER
</span><span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.222"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( ORDER, <span class="string">'B'</span> ) ) THEN
         IORDER = 2
      ELSE IF( <a name="LSAME.224"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( ORDER, <span class="string">'E'</span> ) ) THEN
         IORDER = 1
      ELSE
         IORDER = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for Errors
</span><span class="comment">*</span><span class="comment">
</span>      IF( IRANGE.LE.0 ) THEN
         INFO = -1
      ELSE IF( IORDER.LE.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( IRANGE.EQ.2 ) THEN
         IF( VL.GE.VU ) INFO = -5
      ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) )
     $          THEN
         INFO = -6
      ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) )
     $          THEN
         INFO = -7
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.249"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SSTEBZ.249"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize error flags
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      NCNVRG = .FALSE.
      TOOFEW = .FALSE.
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Simplifications:
</span><span class="comment">*</span><span class="comment">
</span>      IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
     $   IRANGE = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants
</span><span class="comment">*</span><span class="comment">     NB is the minimum vector length for vector bisection, or 0
</span><span class="comment">*</span><span class="comment">     if only scalar is to be done.
</span><span class="comment">*</span><span class="comment">
</span>      SAFEMN = <a name="SLAMCH.274"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> )
      ULP = <a name="SLAMCH.275"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      RTOLI = ULP*RELFAC
      NB = <a name="ILAENV.277"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="SSTEBZ.277"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a>'</span>, <span class="string">' '</span>, N, -1, -1, -1 )
      IF( NB.LE.1 )
     $   NB = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Special Case when N=1
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.1 ) THEN
         NSPLIT = 1
         ISPLIT( 1 ) = 1
         IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN
            M = 0
         ELSE
            W( 1 ) = D( 1 )
            IBLOCK( 1 ) = 1
            M = 1
         END IF
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute Splitting Points
</span><span class="comment">*</span><span class="comment">
</span>      NSPLIT = 1
      WORK( N ) = ZERO
      PIVMIN = ONE
<span class="comment">*</span><span class="comment">
</span>CDIR$ NOVECTOR
      DO 10 J = 2, N
         TMP1 = E( J-1 )**2
         IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
            ISPLIT( NSPLIT ) = J - 1
            NSPLIT = NSPLIT + 1
            WORK( J-1 ) = ZERO
         ELSE
            WORK( J-1 ) = TMP1
            PIVMIN = MAX( PIVMIN, TMP1 )
         END IF
   10 CONTINUE
      ISPLIT( NSPLIT ) = N
      PIVMIN = PIVMIN*SAFEMN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute Interval and ATOLI
</span><span class="comment">*</span><span class="comment">
</span>      IF( IRANGE.EQ.3 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        RANGE='I': Compute the interval containing eigenvalues
</span><span class="comment">*</span><span class="comment">                   IL through IU.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute Gershgorin interval for entire (split) matrix
</span><span class="comment">*</span><span class="comment">        and use it as the initial interval
</span><span class="comment">*</span><span class="comment">
</span>         GU = D( 1 )
         GL = D( 1 )
         TMP1 = ZERO
<span class="comment">*</span><span class="comment">
</span>         DO 20 J = 1, N - 1
            TMP2 = SQRT( WORK( J ) )
            GU = MAX( GU, D( J )+TMP1+TMP2 )
            GL = MIN( GL, D( J )-TMP1-TMP2 )
            TMP1 = TMP2
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         GU = MAX( GU, D( N )+TMP1 )
         GL = MIN( GL, D( N )-TMP1 )
         TNORM = MAX( ABS( GL ), ABS( GU ) )
         GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
         GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute Iteration parameters
</span><span class="comment">*</span><span class="comment">
</span>         ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
     $           LOG( TWO ) ) + 2
         IF( ABSTOL.LE.ZERO ) THEN
            ATOLI = ULP*TNORM
         ELSE
            ATOLI = ABSTOL
         END IF
<span class="comment">*</span><span class="comment">
</span>         WORK( N+1 ) = GL
         WORK( N+2 ) = GL
         WORK( N+3 ) = GU
         WORK( N+4 ) = GU
         WORK( N+5 ) = GL
         WORK( N+6 ) = GU
         IWORK( 1 ) = -1
         IWORK( 2 ) = -1
         IWORK( 3 ) = N + 1
         IWORK( 4 ) = N + 1
         IWORK( 5 ) = IL - 1
         IWORK( 6 ) = IU
<span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLAEBZ.367"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
     $                WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
     $                IWORK, W, IBLOCK, IINFO )
<span class="comment">*</span><span class="comment">
</span>         IF( IWORK( 6 ).EQ.IU ) THEN
            WL = WORK( N+1 )
            WLU = WORK( N+3 )
            NWL = IWORK( 1 )
            WU = WORK( N+4 )
            WUL = WORK( N+2 )
            NWU = IWORK( 4 )
         ELSE
            WL = WORK( N+2 )
            WLU = WORK( N+4 )
            NWL = IWORK( 2 )
            WU = WORK( N+3 )
            WUL = WORK( N+1 )
            NWU = IWORK( 3 )
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
            INFO = 4
            RETURN
         END IF
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        RANGE='A' or 'V' -- Set ATOLI
</span><span class="comment">*</span><span class="comment">
</span>         TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
     $           ABS( D( N ) )+ABS( E( N-1 ) ) )
<span class="comment">*</span><span class="comment">
</span>         DO 30 J = 2, N - 1
            TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+
     $              ABS( E( J ) ) )
   30    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         IF( ABSTOL.LE.ZERO ) THEN
            ATOLI = ULP*TNORM
         ELSE
            ATOLI = ABSTOL
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( IRANGE.EQ.2 ) THEN
            WL = VL
            WU = VU
         ELSE
            WL = ZERO
            WU = ZERO
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
</span><span class="comment">*</span><span class="comment">     NWL accumulates the number of eigenvalues .le. WL,
</span><span class="comment">*</span><span class="comment">     NWU accumulates the number of eigenvalues .le. WU
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      IEND = 0
      INFO = 0
      NWL = 0
      NWU = 0
<span class="comment">*</span><span class="comment">
</span>      DO 70 JB = 1, NSPLIT
         IOFF = IEND
         IBEGIN = IOFF + 1
         IEND = ISPLIT( JB )
         IN = IEND - IOFF
<span class="comment">*</span><span class="comment">
</span>         IF( IN.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Special Case -- IN=1
</span><span class="comment">*</span><span class="comment">
</span>            IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN )
     $         NWL = NWL + 1
            IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN )
     $         NWU = NWU + 1
            IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE.
     $          D( IBEGIN )-PIVMIN ) ) THEN
               M = M + 1
               W( M ) = D( IBEGIN )
               IBLOCK( M ) = JB
            END IF
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           General Case -- IN &gt; 1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute Gershgorin Interval
</span><span class="comment">*</span><span class="comment">           and use it as the initial interval
</span><span class="comment">*</span><span class="comment">
</span>            GU = D( IBEGIN )
            GL = D( IBEGIN )
            TMP1 = ZERO
<span class="comment">*</span><span class="comment">
</span>            DO 40 J = IBEGIN, IEND - 1
               TMP2 = ABS( E( J ) )
               GU = MAX( GU, D( J )+TMP1+TMP2 )
               GL = MIN( GL, D( J )-TMP1-TMP2 )
               TMP1 = TMP2
   40       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            GU = MAX( GU, D( IEND )+TMP1 )
            GL = MIN( GL, D( IEND )-TMP1 )
            BNORM = MAX( ABS( GL ), ABS( GU ) )
            GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN
            GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute ATOLI for the current submatrix
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABSTOL.LE.ZERO ) THEN
               ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) )
            ELSE
               ATOLI = ABSTOL
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( IRANGE.GT.1 ) THEN
               IF( GU.LT.WL ) THEN
                  NWL = NWL + IN
                  NWU = NWU + IN
                  GO TO 70
               END IF
               GL = MAX( GL, WL )
               GU = MIN( GU, WU )
               IF( GL.GE.GU )
     $            GO TO 70
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Set Up Initial Interval
</span><span class="comment">*</span><span class="comment">
</span>            WORK( N+1 ) = GL
            WORK( N+IN+1 ) = GU
            CALL <a name="SLAEBZ.496"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
     $                   D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
     $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
     $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
<span class="comment">*</span><span class="comment">
</span>            NWL = NWL + IWORK( 1 )
            NWU = NWU + IWORK( IN+1 )
            IWOFF = M - IWORK( 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute Eigenvalues
</span><span class="comment">*</span><span class="comment">
</span>            ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
     $              LOG( TWO ) ) + 2
            CALL <a name="SLAEBZ.509"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
     $                   D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
     $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
     $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy Eigenvalues Into W and IBLOCK
</span><span class="comment">*</span><span class="comment">           Use -JB for block number for unconverged eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>            DO 60 J = 1, IOUT
               TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Flag non-convergence.
</span><span class="comment">*</span><span class="comment">
</span>               IF( J.GT.IOUT-IINFO ) THEN
                  NCNVRG = .TRUE.
                  IB = -JB
               ELSE
                  IB = JB
               END IF
               DO 50 JE = IWORK( J ) + 1 + IWOFF,
     $                 IWORK( J+IN ) + IWOFF
                  W( JE ) = TMP1
                  IBLOCK( JE ) = IB
   50          CONTINUE
   60       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            M = M + IM
         END IF
   70 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
</span><span class="comment">*</span><span class="comment">     If NWL+1 &lt; IL or NWU &gt; IU, discard extra eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>      IF( IRANGE.EQ.3 ) THEN
         IM = 0
         IDISCL = IL - 1 - NWL
         IDISCU = NWU - IU
<span class="comment">*</span><span class="comment">
</span>         IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
            DO 80 JE = 1, M
               IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
                  IDISCL = IDISCL - 1
               ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
                  IDISCU = IDISCU - 1
               ELSE
                  IM = IM + 1
                  W( IM ) = W( JE )
                  IBLOCK( IM ) = IBLOCK( JE )
               END IF
   80       CONTINUE
            M = IM
         END IF
         IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Code to deal with effects of bad arithmetic:
</span><span class="comment">*</span><span class="comment">           Some low eigenvalues to be discarded are not in (WL,WLU],
</span><span class="comment">*</span><span class="comment">           or high eigenvalues to be discarded are not in (WUL,WU]
</span><span class="comment">*</span><span class="comment">           so just kill off the smallest IDISCL/largest IDISCU
</span><span class="comment">*</span><span class="comment">           eigenvalues, by simply finding the smallest/largest
</span><span class="comment">*</span><span class="comment">           eigenvalue(s).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           (If N(w) is monotone non-decreasing, this should never
</span><span class="comment">*</span><span class="comment">               happen.)
</span><span class="comment">*</span><span class="comment">
</span>            IF( IDISCL.GT.0 ) THEN
               WKILL = WU
               DO 100 JDISC = 1, IDISCL
                  IW = 0
                  DO 90 JE = 1, M
                     IF( IBLOCK( JE ).NE.0 .AND.
     $                   ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
                        IW = JE
                        WKILL = W( JE )
                     END IF
   90             CONTINUE
                  IBLOCK( IW ) = 0
  100          CONTINUE
            END IF
            IF( IDISCU.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span>               WKILL = WL
               DO 120 JDISC = 1, IDISCU
                  IW = 0
                  DO 110 JE = 1, M
                     IF( IBLOCK( JE ).NE.0 .AND.
     $                   ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN
                        IW = JE
                        WKILL = W( JE )
                     END IF
  110             CONTINUE
                  IBLOCK( IW ) = 0
  120          CONTINUE
            END IF
            IM = 0
            DO 130 JE = 1, M
               IF( IBLOCK( JE ).NE.0 ) THEN
                  IM = IM + 1
                  W( IM ) = W( JE )
                  IBLOCK( IM ) = IBLOCK( JE )
               END IF
  130       CONTINUE
            M = IM
         END IF
         IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
            TOOFEW = .TRUE.
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If ORDER='B', do nothing -- the eigenvalues are already sorted
</span><span class="comment">*</span><span class="comment">        by block.
</span><span class="comment">*</span><span class="comment">     If ORDER='E', sort the eigenvalues from smallest to largest
</span><span class="comment">*</span><span class="comment">
</span>      IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN
         DO 150 JE = 1, M - 1
            IE = 0
            TMP1 = W( JE )
            DO 140 J = JE + 1, M
               IF( W( J ).LT.TMP1 ) THEN
                  IE = J
                  TMP1 = W( J )
               END IF
  140       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( IE.NE.0 ) THEN
               ITMP1 = IBLOCK( IE )
               W( IE ) = W( JE )
               IBLOCK( IE ) = IBLOCK( JE )
               W( JE ) = TMP1
               IBLOCK( JE ) = ITMP1
            END IF
  150    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( NCNVRG )
     $   INFO = INFO + 1
      IF( TOOFEW )
     $   INFO = INFO + 2
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SSTEBZ.649"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
